///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file src/matrix_math.h
///
/// \brief   header file for functions related to matrix mathematics
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_MATRIX_MATH_HEADER
#define PCPS_MATRIX_MATH_HEADER

namespace pcps {

  template<class S> class MantExp;

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   An exception to be thrown when a matrix cannot be inverted
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////

  class SingularMatrixException : public std::runtime_error {

    public:

      // constructor from string
      SingularMatrixException(const std::string &str) : std::runtime_error(str) {}

  };

  template <class S> void matrix_inverse_lu(const int n,
                                            S & det,
                                            S * const A,
                                            S * const work,
                                            int * const iwork);

  template <class S> pcps::MantExp<S> matrix_determinant_lu(const int n,
                                                            S * const A,
                                                            int * const iwork);

  template <class S> S matrix_determinant_lemma(const int n,
                                                const int m,
                                                const S * const Z,
                                                const S * const U,
                                                const S * const V,
                                                const S * const C,
                                                S * const UZT,
                                                S * const VZ,
                                                S * const Q,
                                                int * const Umask,
                                                int * const Vmask,
                                                int * const Qmask,
                                                S * const work,
                                                int * const iwork,
                                                const bool c_eye);

template <class S> void woodbury_matrix_identity(const int n,
                                                 const int m,
                                                 S * const Z,
                                                 const S * const U,
                                                 const S * const V,
                                                 const S * const C,
                                                 S * const UZT,
                                                 S * const VZ,
                                                 S * const Q,
                                                 int * const Umask,
                                                 int * const Vmask,
                                                 int * const Qmask,
                                                 S * const work,
                                                 int * const iwork);

template <class S> void matrix_math_build_VZ_UZT(const int n,
                                                 const int m,
                                                 const S * const Z,
                                                 const S * const U,
                                                 const S * const V,
                                                 S * const UZT,
                                                 S * const VZ,
                                                 int * const Umask,
                                                 int * const Vmask);

template <class S> void matrix_math_build_Q(const int n,
                                            const int m,
                                            const S * const Z,
                                            const S * const U,
                                            const S * const V,
                                            S * const UZT,
                                            S * const VZ,
                                            S * const Q,
                                            int * const Umask,
                                            int * const Vmask,
                                            int * const Qmask);

} // end namespace pcps

#endif
